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ABSTRACT 

We present the first 3D simulations to include the effects of dark matter annihilation feedback during 
the collapse of primordial mini-halos. We begin our simulations from cosmological initial conditions 
and account for dark matter annihilation in our treatment of the chemical and thermal evolution of 
the gas. The dark matter is modelled using an analytical density profile that responds to changes in 
the peak gas density. We find that the gas can collapse to high densities despite the additional energy 
input from the dark matter. No objects supported purely by dark matter annihilation heating are 
formed in our simulations. However, we find that the dark matter annihilation heating has a large 
effect on the evolution of the gas following the formation of the first protostar. Previous simulations 
without dark matter annihilation found that protostellar discs around Population III stars rapidly 
fragmented, forming multiple protostars that underwent mergers or ejections. When dark matter 
annihilation is included, however, these discs become stable to radii of 1000 AU or more. In the cases 
where fragmentation does occur, it is a wide binary that is formed. 
Subject headings: Dark matter, Population III star formation 



1. INTRODUCTION 

Modern cosmology describes a Universe in accelerated 
expansion, with 73% of its energy density today consist- 
ing of an unknown substance dubbed "dark energy" , and 
the remaining 27% consisting primarily of matter, with a 
very minor contribution from radiation. It is remarkable 
that of the whole matter content, only 17% of it can be 
accounted for by known particles described by the stan- 
dard model of particle physics, whereas the remainder 
is of unknown nature, only (at most) weakly interacting 
with its surroundings. The nature of this dark matter 
(DM) is still unknown, and theories such as supersym- 
metry that extend the standard model and provide can- 
didates for the particle making up the DM are currently 
being tested at the Large Hadron Collider at CERN. 

In spite of the unknown nature of the basic constituents 
of our Universe, its description within the standard cos- 
mological model holds solidly. The growth of structure is 
predicted to take place in a hierarchical fashion. Smaller 
DM and gas structures (halos) are predicted to collapse 
gravitationally at earlier times than more massive ones, 
and evidence for this is indeed observed. 

So far, however, the first generation of stars to have 
formed in the Universe, after the long period of cos- 
mic silence and darkness known as the cosmic "dark 
ages" has not been observed. The formation of stars 
from this first generation - Population III stars - is har- 
boured in the bosom of the very first DM and gas halos 
to undergo gravitational collapse: structures of masses 
10 5 < M/M Q < 10 7 which become nonlinear at redshifts 
30 < z < 20. These stars are too distant to be observed 
directly with any current telescope, and so their nature 
must be probed theoretically. Originally, it was thought 
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that i solated, extremely m assive stars (|Abel et al.ll2000L 
120021 : lYoshida et all 120081 ) would be formed from such 
primordial halos as they have high temperatures and 
only simple molecules such as H2 and HP to provide th e 
necessary cooling (e.g. [Glover 2005; Omukai et al.ll2005l) . 
However, recent work has shown even in these conditions 
small stellar systems may be formed due to fragmenta- 
tion occ urring in protostellar discs around Population 
III st ars dClark et al.ll2q08t :lTurk et al.|[2009t IStacv et al.l 



[MaiClark et a.1.11201 lallbt i Greif et al. 1 120111120121) . Such 



systems are predicted to have masses ranging from a few 
tenths to a few tens of a solar mass, with a m ass spectrum 
that is presumably flat (jDopcke et al.ll20T2T ). 

Interestingly enough, the story of the first stars could 
be intimately related to that of the DM. One popular 
class of candidates for the DM particle are the weakly 
interacting massive particles (WIMPs). These are the 
lightest supersymmetric partners in supersymmetry the- 
ories conserving i?-parity, or the lightest mode in extra 
dimension theories under T-parity conservation. They 
are stable and they "naturally" com ply with all the phe - 
nomenological requirements of DM dTaoso et al. 2008a) . 
Although we refer the reader to iBertone et alj (|2005l ) for 
a detailed review of WIMPs, we note here that they have 
two very important properties: they are self-annihilating 
Majorana particles (i.e. they are their own anti-particles) 
and they interact only weakly with baryons. 

While the effects of WIMPs on st ellar formation i n 
the local Un iv erse are suppressed (|Ascasibar.l 2007), 
iSpolvar et alj (|2008[ ) noted that the annihilation of 
WIMPs accumulated by gas-driven gravitational drag 
into the center of a collapsing halo would strongly heat 
the gas at the center of the halo. In the absence of metals 
- i.e. in a halo forming Population III stars - the energy 
input could potentially be so large as to overcome the 
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cooling provided by H2 molecules. The authors specu- 
lated that this could halt the collapse of the gas, and 
form what they dubbed a "dark star": a pseudo-star 
with a size of an AU or more powered by DMA rather 
than by nuclear fusion. It was also soon realised that if a 
canonical protostar were to be formed, it could capture 
DM through weak scattering between the WIMPs and 
the baryons constituting the star, and that this could in 
principle affect the main seque nce (MS) evolution of the 
newly formed star (|Ioccoll2008l ). 

Several studie s have since addres sed t hese processes 
in mo re detail. iFreese et al.l (j2008l ) and iSpolvar et al.l 
( 2009) proposed that the accretion of infalling gas could 
be sustained by DMA indefinitely up to the formation o f 
objects of masses of 10 5 M Q . However. Hocco et al.l (2008) 
found that such a phase could not last for longer than 10 4 
yr before rea c hing t he Hayashi track, and the findings in 
iHirano et all (f201lh were intermediate between these ex- 
tremes. The results of those studying the capture of DM 
onto a hydrostatic object mediated via weak interaction, 
are more homogenous. It is found that weakly captured 
DM can in principle indefinitely power the star as it is 
sustained not by nuclear reactions bu t by DMA lead- 
ing to a prolonga t ion of the star's life (jlocco et alj|2008t 
lYoon et al.l 120081: iTaoso et all l2008b[ ). It is not clear, 
however, how much DM can be consumed in the pre-MS 
process, and therefore how much is available to be cap- 
tured during the main sequence ([Sivertsson fc Gondolol 
|20H . Indeed, some of the more extreme models of 
dark star formation create problems in standard reion- 
isation scenarios due to their predicti ons of very high 
stellar mass es or very long lifetimes (jSchleicher et al.1 
120081 I2009al ). However, in some cases dark stars may 
delay rather tha n accelerate the reionisation process 
(jScott et al.ll201lD . Direct obse rvation seems cha ll engin g 
for most models of dark stars, iZackrisson et al.l (|2010f ). 
whereas the most extreme objects can already be ruled 
out w ith the use of HST observations, (jZackrisson et al.l 
l2010h . 

Regardless of the efficiencies of these processes, DMA 
takes place throughout the halo. The effects of DMA 
will be felt most keenly when the halo has collapsed 
to high central gas densities, but subtle effects due to 
the influence of annihilation upon the gas chemistry can 
take place ev en earlier in the h alo's life. This has been 
addressed in iRipamonti et all (Hp 10), who used a one- 
dimensional code to self-consistently solve for the chemi- 
cal, thermal and dynamical evolution of the cloud in the 
presence of DMA. Intriguingly, the authors found that 
DMA had little effect on the gas prior to its collapse to 
near-stellar densities. However, the scope of these re- 
sult s was limited as t he ma ximum resolvable density in 
the IRipamonti et al.l (|2010f ) study was reached before a 
hydrostatic object was formed. In complementary work, 
iStacv et"aLl (|2012h addressed the phases after the forma- 
tion of a hydrostatic object, adopting an ad-hoc prescrip- 
tion for the earlier phase. These authors found that if 
multiple protostars were present, dynamical interactions 
would displace them from the DM density peak, thereby 
removing them from the region where DMA would have 
the most influence. 

In this paper we bridge the gap before and after star 
formation in one continuous three-dimensional simula- 



tion that fully accounts for the effects of DMA, primor- 
dial gas chemistry and time-dependent heating and cool- 
ing. This allows us to determine at which stage, if any, 
DMA plays a role in the formation of a Population III 
star. Our paper is structured as follows: Section [2] out- 
lines the method, Section |3] describes the initial condi- 
tions, Sections H] and [S] present our results, Section [S] dis- 
cusses the results, and finally Section [7] summarises our 
conclusions. 

2. METHOD 
2.1. Basic approach 

We perform the calculations for this paper using 
the smo othed particle h ydrodynamics (SPH) code GAD- 
GET 2 (jSpringell [2005) . We have substantially modi- 
fied this code to include a fully time-dependent chem- 
ical network, details of which can be found in the 
Appendix of iClark et all (1201 laf). Our trea tment in- 
cludes: H2 line cooling (Gl over fc Abel I2008D . which is 
treated in the optically thick regime using the Sobole v 
approximation (jYoshida et al.ll2006t IClark et al.ll2011al ): 
collision-induced emissio n from very high density H2 
(jRipamonti fe A bel 2004); cooling due to the collisional 
ionisation and recombina tion of hydrogen and helium 
(Gl over fc Mac Lowll2007l) ; as well as heating and cooling 
arising due to changes in the molecular fraction, shock 
heating, c ompressiona l heat ing and cooling due to rar- 
efactions. iTurk et all (|2011l ) showed that the choice of 
the H2 three-body formation rate coefficient influences 
the resulting dynamics of the gas within the halo. In 
this work we use the three-body H2 formation rate of 
IGloverl (|2008l ) which is intermediate within the range of 
the published rates. 

To treat gas which collapses to scales smaller than the 
resolution limit of our simulations, we use a s ink particle 
appro ach. We create sinks using the standard lBate et all 
( 1995) algorithm. Our parti cular implemen t ation of this 
is the same as that used in Uappsen et al.l ([2005) . The 
minimum number density for sink particle creation is 
10 16 cm -3 , but it is important to note that candidate 
regions must also satisfy a series of tests to ensure that 
they are unambiguously collapsing. The ratio of ther- 
mal energy and rotational energy to the magnitude of 
the gravitational potential energy (a and /3 respectively) 
for the particles that will be turned into the sink must 
satisfy the conditions 

a<\, (1) 
a + (3<l. (2) 

In addition, the total energy of the particles must be 
negative, and the divergence of their accelerations must 
be zero in order to prevent structures that are likely to 
bounce or be tidally disrupted being turned into sinks. 
Consequently sink particle formation often occurs above 
this density threshold. 

We set the outer accretion radius of the sink particles, 
?"out! to 6 AU and the inner accretion radius, r m , to 4 AU. 
SPH particles that come within a distance r out > r > r ln 
of a sink are accreted only if they are gravitationally 
bound to that sink. SPH particles that come closer than 
n n are always accreted. Once sinks have formed, we ac- 
count for the accretion luminosity produced by the ongo- 
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ing accret ion of gas on t o the sinks using the scheme de- 
scribed in Sm ith et al.l ([201 ll ). which assumes that all of 
the protostars that form are normal Population III pro- 
tostars, and not dark stars. We would expect the heating 
produced by accretion onto dark stars to be smaller, ow- 
ing to their larger radii, and so this procedure gives us 
an upper limit on the effect of the accretion luminosity. 

2.2. Treatment of dark matter annihilation 

Most current simulations have insufficient resolution to 
follow the DM contraction within a primordial halo with 
the same resolution as tha t used for the baryons. For 
example, lAbel et al.l (|2002ft had a DM mass resolution 
of 1.1 M Q which meant they could not determine the 
DM profile at radii smaller than around 0.05 pc. In- 
ste ad, the method o f adiabatic contraction developed 
by iBlumenthal et al.l (119861) has been used by vari- 
ous a uthors (e.g. iSpolvar et al.l 120081 : iRipamonti et al.1 
2010) to calculate how the DM distribution will re- 
spond to changes in the gravitational potential of the gas. 
ISpolvar et al.l (|2008f ) a pplied this method to t he collaps- 
ing ha lo simulations of lAbel et al.l ([20021 ) and lGao et al.l 
(2007) and found that the DM density at the outer edge 
of the baryonic core after adiabatic contraction was well- 
fit by the expression 

p xc 5n° p 81 GeV cm" 3 , (3) 

where p xc is the DM density at the edge of the core and 
n p is the number density of hydrogen nuclei in the coreQ 
Given the considerable technical challenges involved if 
one attempts to follow the evolution of the DM density 
profile directly, we do not attempt to do so. Instead, 
we parameterise its effects, with the help of th e adia- 
batic contraction results of ISpolvar et al.l (|2008l ). The 
DM density profile is calculated within the code as fol- 
lows. At each time-step, the point of maximum gravi- 
tational potential energy in the halo is found, which we 
assume to be the point at which the DM density has its 
maximum value. Next, the maximum baryon density of 
the halo is found and the DM density at the outer edge 
of the core calculated from it using Equation [3] Its ra- 
dial dependence is then described by a two-part power 
law fit, using power-law s l opes drawn from the calcula- 
tions of IRipamonti et all ((2010). In the outer halo, at 
distances r > r c from the centre of the density profile, 
we have 

Px = p0 1 r~c" ) GeVcm ~ 3 < ( 4 ) 

where r is the radial distance from the centre and po — 
5 x 10 4 GeV cm -3 is a normalising constant that is equal 
to the DM density at r = 1 pc. The size of the central 
DM core, r c , is determined by finding the radius at which 
the DM densities given by Equations [3] and [4] are equal, 
which yields 

/ n \ -0.81/1.8 

r - = 16 ' 7 (io^) au - < 5 > 

1 To convert from the DM density in energy units to the value 
in the more familiar mass density units, one simply multiplies by a 
factor of c 2 . Moreover, it is also possible to express the gas density 
in energy units; for a gas with primordial composition, we have 
p g ~ 1.24 n GeV cm -3 , where n is the number density of hydrogen 
nuclei. 



Finally, the DM density profile within the core, at radii 
r < r c , is given by 

( r \ ~°- 5 

Px=Pxc\— J ■ (6) 

It is unclear how efficient adiabatic contraction will be 
within the central regions of the DM halo, given that 
gravitational collapse of the gas in these regions occurs 
rapidly. Our adoption of a peaked DM profile within the 
core may therefore be an overestimate; see, for instance, 
the discussion of the val idity of the adiabatic co ntraction 
approximation given in IRipamonti et all (|2010| ) . 

Given the DM density profile, the contribution of DMA 
to t he halo energy budge t is computed as follows. We fol- 
low ISpolvar etaLl ([2008D and adopt the standard DMA 
cross section of (av) = 3 x 10~ 26 cm 3 s _1 . For the DM 
particle mass, m x , we note that astrop hysical constraints 
require that m x > 10 GeV (see e.g. lAckermann et al.l 
l201lHGalli et al.ll2011ft . Additionally the desire to avoid 
fine-tuning while still explaining the gauge hierarchy 
problem argues for an upper limit m x ~ a few TeV. We 
therefore adopt a DM particle mass of m x = 100 GeV 
for our fiducial case, close to the middle of this range 
of values, but we also explore the effect of varying m x . 
Following Val des fc Ferr ara (20Q8|), we assume that one 
third of the energy from DMA is immediately carried 
away by neutrinos and that the remaining energy is split 
between direct heating of the gas, ionisation and disso- 
ciation of the constituents of the gas, and the excitation 
of H, He and H2. The fraction that is deposited as heat 
is given by 

f h = 1 - 0.875(1 - z^ 4052 ), (7) 
and the fraction that goes into ionisation is given by 

fi = 0.384(1 _ a; 0.542 ) 1.1952 j (8) 

where x e = n c /n is the fractional abundance of electrons. 
The remaining energy is radiated away in electronic tran- 
sitions from H, He and H2 and plays no direct role in the 
evolution of the gas. The total power per unit volume 
injected by DMA is 

W ann = 2m x N aim , (9) 

where -/V ann is the number of DMAs per unit volume per 
unit time. This is given by 

N,nn=\n 2 x (<7v), (10) 

where n x = p x /m x is the number density of DM parti- 
cles. The resulting energy injection rates per unit volume 
for heating (Qh) and ionisation (Qi) are 

Qh = 2/ a A/ann(l - 6^)/^, (11) 

Qi = 2/ a A/ ann (l - e- T *)f imx , (12) 

where f a = 2/3 is the fraction of the total energy that 
affects the gas (one third of the energy escapes in the 
form of neutrinos) . 

We estimate the optical depth of the gas to the an- 
nihilation products, t x , by assuming that the baryonic 
density profile is fiat within the DM core and declines 
as a power-law p cx r~ 2 at larger radii. For a profile of 
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TABLE 1 

Overview of simulations 



Simulation 


Annihilation 


DM particle mass [GeV] 


Hl-ref 


no 





Hl-lm 


yes 


10 


HI 


yes 


100 


Hl-hm 


yes 


1000 


H2-ref 


no 





H2 


yes 


100 



this form, the column density of gas, measured radially 
outwards from a point r, is given at r < r c by 



S(r) = p(r) 
and at r > r c by 



n 



E(r) = p(r) [ r - — ) , 



(13) 



(14) 



where r n is the virial radius of the halo, which we take 
to mark the edge of t he gas distribution. If we follow 
iRipamonti et al.l(|2010ft and adopt a constant gas opacity 
of k = 0.01 cm 2 g _1 , and also restrict our attention to 
scales small compared to r/j, then we can write the optical 
depth as 



kS(t) 



np(r)(2r c 
np(r)r 



r) r < r c 
r > r c 



(15) 



In practice, this expression for t x is an underestimate, 
for a couple of reasons. First, it assumes that the opti- 
cal depth in any direction from point r is the same as 
the value that we measure along a ray passing radially 
outwards, when in reality the optical depths in other di- 
rections will be somewhat larger. Second, the baryonic 
density profiles that we recover in our simulations do not 
completely match up with the profile we assumed to de- 
rive t x . Our assumed profiles and the real profiles both 
have power-law slopes p oc r~ 2 in their outer regions, 
but the real profiles do not have the completely flat core 
that we assume within r c . Our estimate of t x is therefore 
too small at distances r >C r c , meaning that our derived 
heating and ionisation rates may also be too small. In 
practice, we do not expect this to be a major source of 
error, as in the regime where DMA becomes important, 
gas with r r c will typical have t x > 1, meaning that 
the heating and ionisation rates are insensitive to the 
precise value of t x . 

To convert from Qi, the total energy per unit time per 
unit volume that goes into ionisations, to the ionisation 
rates of the individual chemical speci es present in the 
gas, w e follow the procedure outlined in IRipamonti et al.l 
(2007). They identified seven main reactions that can be 
caused by DMA: the ionisation of H, He, He + and D, and 
the dissociation of H2, HD and H^. Between them, these 
seven reactions account for the bulk of the energy lost 
in the form of ionisatio ns or dissociations, and we follow 
IRipamonti et al.l (|2007j) and split the energy available for 
ionisation between these seven species according to their 
relative abundances. 

3. INITIAL CONDITIONS 

We take our initial conditions fr om the cosmolog i- 
cal simulations and resimulations by IGreif et al.l (|201lD . 



These simula tions made use of the novel moving mesh 
code AREPO (|Springelll2010D to fully resolve the forma- 
tion of five minihalos from cosmological initial conditions. 
Cells were refined during the evolution to ensure that the 
Jeans length was always resolved by at least 128 mesh 
points, up until the point at which the maximum num- 
ber density in the collapsing gas reached n = 10 9 cm~ 3 . 
All of the' halos modelled by IGreif et all (|2011[) formed 
multiple protostars with a range of masses. 

For this work w e cut out the central two parsecs of the 
IGreif et "all (|2011| ) simulations and continue their evolu- 
tion using our modified version of GADGET 2 with DMA, 
implemented as discussed in the previous section. The 
halos were selected at the point where their central gas 
number density reached n = 10 6 cm -3 for the first time. 
We selected this point to begin our resimulation because 
preliminary modelling using simplified initial conditions 
showed that this was the point at which indirect feedback 
from DM A- induced ionisation first becomes important. 

Each mesh point in AREPO is converted to an SPH par- 
ticle with the same properties (mass, momentum, etc.) 
as the original mesh cell associated with that mesh point. 
As the network for primordial chemistry that is imple- 
mented in AREPO was derived from the network that 
we use in GADGET, both codes evolve the same set of 
chemical species and it is straightforward to transfer the 
chemical abundances from one code to the other. The 
formation of the first protostar occurs in the central re- 
gion of the halos where the SPH particle masses are 
10~ 4 Mm, giving us a mas s resolution of around 10~ 2 
M Q (jBate fc Burkertlll997D . For a test of the effective- 
ness of using AREP O initial conditions for gadget 2 see 
ISmith et all (|20T1. 



IGreif et al.1 (2011) simulated the evolution of five dif- 
ferent DM halos. In this work, however, we fo cus on only 
two o f these five: Halo 1 (in the notation of IGreif et al.l 
120 111) , which in the original calculation rapidly fragments 
into a multiple system, and Halo 2, which undergoes a 
phase of HD-dominated co oling, which ultima tely leads 
to less fragmentation (c.f. iClark et al.1 120 llal ). Halo 1 
has a mass of 1810 M Q within its central 2 pc and we 
denote it in our study as HI. Halo 2 has a mass of 1240 
M Q within its central 2 pc, and will be referred to as H2. 
Halo 1 has only ~ 6.9 x 10 5 SPH particles and Halo 2 
~ 6.3 x 10 5 particles, but due to the on-the-fly refinement 
used in AREPO the particle mass in the central regions of 
the halo is only 10 -4 M Q as mentioned above. We did 
not include any traditional dark matter particles, but in- 
stead treated the dark matter analytically as described 
in the previous section. 

For both halos, we run one simulation in which the 
effects of DMA are not included (Hl-ref, H2-ref) and a 
second which assumes a DM particle mass of 100 GeV 
(HI, H2). For Halo 1, we also run two simulations with 
different values of m x : one in which we set m x = 10 GeV 
(Hl-lm) and a second in which we set m x = 1000 GeV 
(Hl-hm). As the power produced per unit volume by 
DMA scales as m x n x oc m x p^./m x oc p x /m x , these two 
runs correspond to cases in which the energy input rate is 
increased or decreased by an order of magnitude, respec- 
tively. As astrophysical constraints st rongly disfavour 
DM masses smaller than m x = 10 GeV, ([S chleiche r et al.l 
l2009bHAckermann et al.ll20ll IGalli et al.ll2011h . our HI- 
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TABLE 2 
Time taken to form first sink 



Simulation 


Time (10 s yr) 


Hl-ref 


3.57 


HI 


2.38 


Hl-lm 


2.70 


Hl-hm 


2.57 


H2-ref 


3.59 


H2 


2.36 



lm model should give an indication of the largest effect 
that we can reasonably expect to obtain from DMA. De- 
tails of our simulations are summarised in Table [T] 

4. COLLAPSE TO NEAR STELLAR DENSITIES 
4.1. Densities and temperatures 

In Figure [TJ we show the density profiles of gas and 
DM for simulations HI, Hl-lm, and H2, plus the two ref- 
erence models. The simulations are compared when the 
hydrogen nuclei number density at the centre of the halo 
first reaches 5 x 10 14 cm -3 , which corresponds to a dark 
matter core radius of r c ~ 8 AU and a time of roughly 
6 years before the formation of the first protostar. It is 
immediately apparent that the gas in each halo is able 
to collapse to high densities, regardless of the strength 
of the DMA feedback. This is true even in our maxi- 
mal feedback model, Hl-lm, where the DM particle mass 
was only 10 GeV. For reference, in lSpolvar et al.l (|2008f ) 
it was predicted that collapse would stop at densities of 
10 13 cm -3 for a DM mass of 100 GeV, and at densi- 
ties of 10 9 cm" 3 for a DM mass of 10 GeV. We find 
no evidence for this in our simulations. On the other 
hand, o ur results are in good a greement with the ID re- 
sults of lRipamonti et al.l (|2010f ) . who found that the gas 
could collapse to densities of 10 13 -10 14 cm -3 in all of 
their models. We also found collapse up to high densi- 
ties in run Hl-hm, performed with a DM particle mass 
of 1000 GeV, but as the effects of DMA in this case are 
much smaller than in our fiducial case, we do not discuss 
this run further in this section. 

The main effect that DMA has on the density profile of 
the gas is the clear density enhancement at radii of 100- 
1000 AU. We also find that in the runs with DMA, the 
time taken for the gas to collapse to protostellar densities 
is shorter than in the reference runs (see Table ^ . Both 
of these effects can be understood by an analysis of the 
temperature structure of the gas. In Figure [2] we show 
how the temperature of the gas varies as a function of 
the gas number density In the outer regions of the halo, 
the gas is cooler than in the reference case. At higher 
densities, however, there is a sharp rise in the gas tem- 
perature, taking it above the value in the reference case. 
Comparison of Figures [T] and [2] shows that the density 
at which this temperature increase occurs is the same as 
the density at which we first see the bump in the density 
profile. 

In Figure [3J we show how the fractional abundance^ 
of H2 and e _ vary as a function of density within our 
different models. At low densities, DMA-induced ioni- 

2 Note that these are defined here as the ratio of the number 
density of the species of interest to the number density of hydrogen 
nuclei. This means that gas in which all of the hydrogen is in 
molecular form has an H2 fractional abundance xjj 2 = 0.5. 



sation enhances the number of electrons by two or more 
orders of magnitude. This increases the rate of H2 for- 
mation relative to the case without DMA by increasing 
the effectiveness of the following reaction chain: 

H + e-^H-+ 7 , (16) 
H-+H^H 2 +e-. (17) 

As H2 is the main coolant of the gas at these densi- 
ties, the enhanced H2 abundance leads to more efficient 
cooling. We therefore find, counter-intuitively, that this 
lower density gas is cooler in the case with DMA than 
in the case without DMA. This is in agreement with 
the re sults of the ID simulations of Ripamont i et al.l 
(2010), and provides an explanation for the shorter col- 
lapse times of the halos in which DMA is occurring, as 
for the majority of its lifetime, the gas has less thermal 
support. 

Once the gas reaches a density of around 10 12 cm -3 , 
however, the H2 fraction decreases sharply in the runs 
with DMA. By comparing this Figure with the temper- 
ature distribution shown in Figure [2] we see that this 
sharp decrease occurs once the gas temperature reaches 
T ~ 2000 K. It occurs because at this temperature, col- 
lisional dissociation of H2 by the reactions 



H 2 + H^H + H + H, (18) 

H 2 +H 2 ^H + H + H 2 , (19) 

H 2 +Hc^H + H + He, (20) 

and by the charge transfer reaction 

H 2 + H+ -»• H+ + H, (21) 



becomes effective. The increased temperature leads to 
an increased thermal pressure, which slows the collapse 
of the gas at these densities and leads to a "pile-up" of 
material visible as a pronounced bump in the gas density 
profile. It is also clear from Figure |3J that not all of the 
H2 is destroyed at this density. The H2 that survives in 
the higher density gas enables it to dissipate much of the 
energy introduced into the gas by the DMA, and hence 
allows it to maintain its temperature at close to 2000 K. 
This allows gravitational collapse to continue past the 
point at which DMA heating first outstrips H2 cooling. 
We examine the physics of this in more detail in the next 
section. 

Finally, we note that one must take care in interpret- 
ing the behaviour of the H 2 abundances shown in Fig- 
ure [3J It is all too easy to think of the evolution of xn 2 
with density as a time history of the H2 abundance. In 
this picture, the H2 abundance first increases as the gas 
collapses, then sharply decreases once the number den- 
sity reaches 10 12 cm -3 , before increasing once more at 
higher gas densities. However, this is incorrect, as Fig- 
ure S] demonstrates. The rise in the H2 abundance in the 
densest gas that we see if we look at a single snapshot is 
not due to the reformation of H2 in this gas. Instead, it 
occurs because the amount of H2 that the gas loses once 
it reaches a density of around 10 12 cm' 3 increases as the 
collapse proceeds. The first gas to collapse - the mate- 
rial that is now at densities above 10 14 cm -3 - loses only 
a relatively small fraction of its H2, while the material 
falling in later loses much more of its molecular content. 
In part, this behaviour is due to an increase in the DM 
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Fig. 1. — Density profiles of the gas and the DM at the point at which the hydrogen nuclei number density of the gas at the centre of 
the collapsing core first reached 5 X 10 14 cm -3 . The solid line shows the radially-averaged baryon density and the dotted line shows the 
DM density. The background grey scale shows the baryon density profile before averaging. The dashed line shows the radially-averaged 
baryon density in the reference case without DMA. The main effect of DMA is to create an enhancement in the density profile in the region 
between 100 AU and 1000 AU. 
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Fig. 2. — Temperature of the gas as a function of density, plotted for the time when the hydrogen nuclei number density of the gas at 
the centre of the core first reached 5 X 10 14 cm -3 . The solid line shows the radially-averaged gas temperature at each density and the 
background grey scale shows the gas temperature before averaging. The dashed line shows the radially-averaged gas temperature in the 
reference case without DMA. With DMA, the gas is cooler in the outer, less dense regions of the halo, but much hotter in the dense interior. 
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Fig. 3. — Fractional abundances of H2 and electrons in Halo 1 in the fiducial case (run HI; left panel) and the 10 GeV DM particle mass 
case (run Hl-lm; right panel). The solid lines show the radially-averaged fractional abundances at each density and the background grey 
scale shows the abundances before averaging. The dashed lines show the radially-averaged abundances in the reference case without DMA. 
DMA-induced ionisation increases the electron abundance of the gas at all densities, which in turn leads to an increase in the rate of H2 
formation at low densities. At high densities, there is a drop in the H2 abundance due to DMA-induced dissociation of H2. 
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Fig. 4. — Radially-averaged fractional abundance of H2 as a func- 
tion of density in run HI at three different output times, corre- 
sponding to peak densities of approximately 10 cm -3 (dashed 
line), 10 13 cm -3 (dotted line) and 5 X 10 14 cm -3 (solid line). We 
see that as time goes on, the gas loses an increasingly large fraction 
of its H2 content once it reaches n ~ 10 12 cm -3 . The increase in 
the H2 fraction at higher densities is therefore a consequence of the 
fact that gas which collapsed earlier lost less H2 than gas which 
collapsed later, and does not indicate reformation of H2 within the 
densest gas. 

number density associated with the n = 10 12 cm -3 gas, 
caused by the increasing concentration of the halo. How- 
ever, an additional contribution comes from the fact that 
gas falling in at later times will generally have a higher 
infall velocity, and hence will shock more strongly once 
it reaches this density. 

4.2. Heating and cooling rates 
4.2.1. Fiducial case 

Figure [5] shows the main heating and cooling processes 
acting in run HI shortly before the formation of the first 
protostar, at the time when the central number density 
first reached 5 x 10 14 cm~ 3 . At low densities, comprcs- 
sional heating and H2 formation heating are the two most 
important heat sources, although even at densities as low 
as 10 6 cm~ 3 , DMA heating is beginning to contribute to 
the total heating rate at the level of a few percent or 
more. The cooling of the gas at these densities is dom- 
inated by H2 line cooling. Once the gas density reaches 
n ~ 10 8 cm -3 , compressional heating starts to become 
unimportant, and DMA heating catches up with H2 for- 
mation heating, with both subsequently playing impor- 
tant roles in the thermal balance of the gas. The heating 
produced by these two processes raises the gas temper- 
ature, but this enables H2 line cooling to become more 
effective, and this is able to offset much of the additional 
heat input from the DMA at densities n < 10 12 cm~ 3 . 
Above this density, H2 line cooling becomes relatively 
ineffective, both because the H2 is starting to dissociate 
and also because the H2 emission lines become optically 
thick at high densities. The heating rate due to DMA 
therefore outstrips the H2 line cooling rate, i n agreement 
with the prediction of iSpolvar et al.l (|2008l ). However, 
we see from Figure [5] that this is not the whole story. As 



the temperature rises, H2 collisional dissociation (reac- 
tions [18rt2"0"|) and the destruction of H2 by charge trans- 
fer with H + (reaction [21"]) become increasingly effective. 
These reactions are endothermic, and hence remove ther- 
mal energy from the gas. If H2 subsequently reforms, 
then the associated H2 formation heating will return 
some of this energy to the gas, but if more H2 is destroyed 
than can reform, then the overall effect is to dissipate en- 
ergy. This effect should not be regarded as "cooling" in 
the same sense as H2 line cooling or CIE cooling, as it 
will not bring about a decrease in the gas temperature. 
Instead, it acts more like a thermostat, preventing the 
temperature from increasing significantly until all of the 
H2 has been consumed. 

The importance of this effect can be appreciated if we 
compare the time required to destroy all of the H2 at 
the centre of the halo with the free-fall collapse time of 
the gas. By considering the energy budget of the gas 
in this way we can obtain a rough estimate of whether 
collapse is likely to be halted at densities close to our 
resolution limit. From Figure [5j we see that at a density 
n = 5x 10 14 cm -3 , the DMA heating rate per unit volume 
at the centre of the halo is roughly 5 x 10~ 6 ergs _1 cm~ 3 . 
If we take the H2 fractional abundance in this gas to 
be ih 2 = 0.4, which is a conservative estimate, then 
the total amount of energy per unit volume that can be 
dissipated by dissociating H2 is 

En 2 =4.48eV x x H2 n, (22) 

- U00 (^VTm =3) erg cm" 3 , (23) 

\5 x 10 14 cm 6 J 

where 4.48 eV is the binding energy of a single H2 
molecule. The time required to destroy all of the H2 
is therefore t^ s = 1400/5 x 10 -6 ~3x 10 8 s, or around 
10 years. In comparison, the free-fall time of the gas at 
this density is around 2 years. The energy input from 
the DM therefore cannot destroy the H 2 rapidly enough 
to prevent the gas from collapsing. At even higher den- 
sities, the DMA heating rate will be larger. However, 
the heating rate per unit volume within the core of the 
DM density profile scales with the central gas density as 
n 1,62 , and hence idis oc n~ - 62 . In comparison, the free- 
fall timescale scales as iff oc n~°- 5 . Therefore, a large 
increase in density is required in order to significantly 
alter the ratio of the H2 dissociation timescale to the 
free-fall timescale. 

In reality collapse of the dense gas is likely to occur 
on a timescale that could be a factor of a few longer 
than the free fall timescale. However, even in this case 
the collapse timescale is shorter than the H2 dissociation 
timescale, and will remain so until the density increases 
significantly. It therefore seems unlikely that collapse 
will halt at a density just above our sink creation thresh- 
old, although we cannot say at exactly what point the 
collapse will halt without performing much higher reso- 
lution simulations. 

Figure [5] shows the main heating and cooling processes 
acting in the reference case, run Hl-ref. In the ab- 
sence of DMA heating, the main heating term is com- 
pressional heating over most of the range of densities 
examined, with H2 formation heating becoming impor- 
tant at densities between 10 10 cm -3 and 10 12 cm -3 , and 
at n > 10 14 cm -3 . The main source of cooling at 
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Fig. 5. — Rates of the major heating and cooling processes active in run HI at the time when the central density first reaches n = 
5 X 10 14 cm -3 . DMA heating is an important contributor to the net heating between densities of 10 8 — 10 12 cm -3 . However, over most of 
this range, it is largely balanced by H2 line cooling, and provides insufficient heating to halt the collapse. Above n ~ 10 12 cm -3 , H2 line 
cooling becomes ineffective, and most of the energy introduced into the gas by DMA heating is dissipated by H2 collisional dissociation 
and the destruction of H2 by charge transfer. 
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Fig. 6. — Rates for the major heating and cooling processes acting in simulation Hl-ref at a point just before the formation of the first 
protostar. We show the rates at a later time than in Figure[5]so that the behaviour in the high density regime dominated by H2 dissociation 
cooling is clear. If we compare the results here with those in Figure [5] we see that H2 dissociation cooling only becomes important at 
n ~ 10 14 cm -3 , in contrast to n ~ I0 11 cm -3 in run HI. This is because the gas temperature at these densities is lower in run Hl-ref than 
in run HI, owing to the absence of DMA heating. 
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n < 10 14 cm' 3 is H2 line cooling, while at higher densi- 
ties, H2 collisional dissociation plays an important role 
in regulating the temperature of the gas. This is in rea- 
sonable agreement with the re sults of other models; for 
instance, Yoshida ct al. (200j|) find that H 2 dissociation 
becomes significant at densities of ~ 10 15 cm -3 , at which 
point the gas temperature is roughly 2000 K. The main 
effect of the DMA heating seems to be simply to bring 
the gas to this state at an earlier point in its evolution. 

4.2.2. Maximal Case 

The previous analysis is for the fiducial case where the 
DM particle mass was 100 GeV. However, even in our 
maximal case, where the DM particle mass was 10 GeV, 
we find broadly similar behaviour. In Figure [71 we show 
the rates of the main heating and cooling processes in 
run Hl-lm shortly before the formation of the first pro- 
tostar. In this case, we see that DMA heating becomes 
important at an earlier time and that H2 dissociation be- 
comes the dominant energy dissipation mechanism once 
the gas density reaches a lower density, n ~ 10 10 cm~ 3 , 
than in the fiducial model. This is consistent with the re- 
sults plotted in Figure^ which show that the gas reaches 
T ~ 2000 K slightly earlier in its evolution. However, the 
subsequent behaviour of the gas is very similar in both 
models. H2 dissociation is such an effective thermostat 
that an increase in the heating rate by a factor of ten 
produces only a small increase in the gas temperature. 
Repeating our previous analysis of the H2 dissociation 
timescale, we find that in this case idis ~ iff, meaning 
that the gas will probably lose most of its H2 before it 
collapses to protostellar densities. It is therefore possible 
that in this extreme case, a "dark star" will form, but if 
so, then its size will be smaller than our mi nimum spatial 
resolu tion of a few AU. For comparison, ISpolvar et al.l 
(2008) predict that even in the 100 GeV dark 
star of size ~ 20 AU should form, and expect that reduc- 
ing the DM particle mass will lead to an even larger dark 
star. We find no evidence for such large DMA supported 
structures in our simulations. 

5. SECONDARY FRAGMENTATION 

So far we have only considered th e initial collapse o f 
the halo to stella r densit ies. However. IStacv et alj (2010) 
and IClark et al.l (|2011b[ ) have shown that the protostel- 
lar accretion disc that builds up around the first pro- 
tostar rapidly becomes unstable and fragments. The 
natural implication is that the first stars to form were 
generally part of multiple systems, which has profound 
consequences for our understanding of primordial star 
formation. It raises the possibility of the fir st stars be 



ing e jected while they were still low mass (IGreif et al 



2011 ) or growing in ma ss through mergers ([Smith et al 



2012t IGreif et"aLl I2012D . Evolving close binary systems 
are also a p otential mechanism fo r creating early gamma 
ray bursts (jBromm fc L ocb 2006) . It is therefore impor- 
tant to discover whether primordial accretion discs are 
still unstable in the presence of DMA. 

In order to enable us to follow the evolution of our 
simulated halos beyond the point at which the first pro- 
tostar forms, we use a sink particle treatment, as outlined 
in Section [2~T1 During the evolution of the halo with sink 
particles we fix the location of the DM peak to exactly 



coincide with the first sink formed. This mimics the ef- 
fect of the DM being 'locked in' to the first star to form. 

Figure [8] shows a comparison of the column density of 
the disc in run HI and in the corresponding reference 
case without DMA, run Hl-ref, at a time 500 years after 
the formation of the first protostar. In the reference case, 
the disc has already fragmented, forming four additional 
protostars. In the case with DMA heating, however, the 
disc has not fragmented. We continued run HI until the 
primary mass reached 15 Mq, which occurred at a time 
t = 17, 527 years after the formation of the first proto- 
star, but saw no disc fragmentation during the whole of 
this period. As our models do not include the effects of 
ionisation feedback from massive stars, we stopped our 
simulation at this point. However, we note that studies 
that have investigated the effect of this ionisation feed- 
back find that it rapidly shuts off the supply of material 
to the disc, an d hence will act to suppress fragmentatio n 
at later times ([Stacy et al.l l2012; Ho sokawa et al.ll201lD . 

An indicati on of disc stab ility is given by the Toomre 
Q parameter ([Toomrel [l964). which has the form 



Q 



ttGS 



(24) 



,where c s is the sound speed, k is the epicyclic frequency, 
and E is the surface density of the disc. In the classical 
Toomre analysis, discs are stable if they have Q > 1. 
Such an analysis is technically only valid when the disc 
mass is much smaller than that of the central protostar, 
and the disc is infinitely thin, neither of which are the 
case here. However, Q still proves to be a reliable guide as 
to whether or not the disc is stable, and helps to demon- 
strate why the inclusion of DMA heating has a large 
effect on the disc stability. Figure [9] shows the radially 
averaged properties of the disc surrounding the central 
protostar for Halo 1. We plot three cases, the reference 
case without DMA (Hl-ref), the fiducial case with DMA 
(HI), and the "minimal" case with a high DM particle 
mass (Hl-hm). The comparison between the disc prop- 
erties is made at a time of 217 yr after the first sink 
particle forms in each simulation. For the reference case, 
this corresponds to a time just before disc fragmentation. 
(The model of low DM mass, Hl-lm, shall be discussed 
in more detail later.) 

All of our discs rotate at a similar rate, for instance 
HI has a mean radial velocity of 4.9 kms -1 in the gas 
denser than 1 x 10 10 cm -3 and Hl-ref has a mean radial 
velocity of 4.2 kms _1 in the gas at this density. Conse- 
quently, the most relevant terms are the surface density 
of the disc and the sound speed. Figure [S] shows that the 
surface density of the disc is not substantially changed 
by the DMA. However, the disc temperature is consider- 
ably increased in the DM case compared to the reference 
case. This raises the sound speed and causes the Q value 
to rise above 1 in the DM cases, whereas in the reference 
case it falls below one and the disc fragments. Figure [§] 
also shows the H2 fraction within the disc. In the refer- 
ence case, the disc is fully molecular, but in the DM cases 
the disc is partially dissociated. As the disc evolves the 
H 2 fraction further decreases until all the molecular gas 
is destroyed. After this point the energy can no longer 
be absorbed by the destruction of H2 and the temper- 
ature of the disc will rapidly rise. Thus it becomes in- 
creasingly difficult for the disc to fragment as it evolves, 
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Fig. 7. — As FigureO but for run Hl-lm. In this case, the heating from DMA is important at all densities. 
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Fig. 8. — Column density of a slice through the midplane of the disc in the central regions of runs HI and Hl-ref at a time 500 years 
after the first sink formed in each case. The slice is 1000 AU thick. In the reference case without DMA, four fragments are formed. In the 
DMA case, there is no secondary fragmentation even at later stages of the evolution. 
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Fig. 9. — Properties of the disc 200 years after the first sink forms. The solid line shows the reference case Hl-ref, the dashed line our 
fiducial case HI, and the dotted line shows our minimal case Hl-hm. In the reference case, the disc has become Toomre unstable, but in 
the two cases with DMA, the disc remains stable. Note that the low DM mass case, Hl-lm does fragment at radii of order 1,000 AU later 
in the simulation. 
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TABLE 3 

PaOTOSTELLAR SEPARATION AND PRIMARY MASS AT THE POINT 
WHEN A SECONDARY PROTOSTAR FORMS 



Simulation 


Separation [AU] 


Primary Mass [Mq] 


Hl-ref 


31.7 


1.63 


HI 






Hl-lm 


1629 


11.8 


Hl-hm 






H2-ref 


14.8 


0.68 


H2 


1012 


10.3 



which explains why we see no fragmentation in run HI 
despite running the simulation to a central primary mass 
ofl5M Q . 

Figure [TOl shows a column density projection for a slice 
through the centre of run H2 at three different output 
times. In this second sink forms at large radii 

where the gas is cooler. The fragmentation in run H2 
is substantially different from that seen in the reference 
model. Without DMA, a secondary sink forms within 
155 years at a distance of 14.8 AU from the primary. 
This is followed by additional fragmentation such that 
when the simulation has run for a thousand years, there 
are 7 sinks. In contrast to this, the DM simulation does 
not form a secondary object until t = 6090 years, and the 
fragmentation occurs at a distance of 1012 AU from the 
primary. The two stars then form a wide binary system. 

Table |3] summarises the secondary fragmentation in all 
of our simulations. It is interesting that secondary frag- 
mentation occurs in run H2 but not in run HI for our 
fiducial DM mass case. Figure [T] shows that in run H2, 
a larger bump forms in the density profile than in run 
HI. Similarly, there is also secondary fragmentation in 
Hl-lm at distances of around a thousand AU where there 
was a significant enhancement in the density profile. At 
such large radii the gas is far enough away from the peak 
DM density that DMA cannot dissociate H2 and the gas 
remains fully molecular, allowing it to cool effectively. 
Runs Hl-lm and H2 both had higher DMA levels than 
HI and Hl-hm and consequently a greater density en- 
hancement outside the DM core. (Our normalisation of 
the DM in Halo 2 meant that the ratio of baryonic mat- 
ter to DM was slightly lower in this case.) Consequently 
it can be concluded that DMA suppresses fragmentation 
in the disc close to the primary but it can actually en- 
courage fragmentation in spiral arms at distances of a 
few thousand AU. 

Another striking difference between the cases with and 
without DMA is the mass of the primary at the time that 
the disc first fragments. In the reference case, the pri- 
mary is a low mass object which is evolving adiabatically. 
Its internal structure i s that of an extended ball of gas 
(jOmukai fc Pallall2003f l and it will be prone to interac- 
tions with close neighbours. With DMA, however, the 
primary is around ten solar masses at the point of frag- 
mentation, at which point it undergoes an expansion as it 
redistributes its internal entropy and will shortly there- 
after begin contracting to the main sequence. The later 
evolutionary stage of the primary and greater distance 
from the secondary in the DMA case mean it is unlikely 
that the stars will influence each other significantly dur- 
ing their formation. 

The left and middle panels of Figure [10] are centred on 
the original position at which the sink particle is orig- 



inally formed. It can be seen that there is some drift 
in the position of the sink after the first few thousand 
years. As our DM halo is analytic and cannot dynami- 
cally respond to the baryons we assume that the DM is 
locked into the first protostar and follows this centre as 
its position migrates. Similarly our idealised DM halo is 
oblivious to the formation of spiral arms in the baryonic 
component of the halo. It is conceivable that both the 
drift of the sink particle and the transfer of energy from 
the spiral arms to the DM distribution could heat the 
central DM distribution and somewhat decrease the DM 
central density. Studies with a highly resolved live DM 
halo and high resolution baryonic component would be 
needed to fully address this issue. 

6. DISCUSSION 

6.1. Do dark stars form? 

In their original study, iSpolvar et al.l (|2008[) made the 
simple assumption that the gravitational collapse of the 
gas would come to a halt as soon as the heating rate 
produced by DMA exceeded the H2 cooling rate. As a 
result, they predicted the formation of large (~ 20 AU or 
more) protostellar objects supported by DMA heating, 
the so-called "dark stars". However , their assumption 
is incorrect. iRipamonti et al.l (|2010l) first showed, us- 
ing ID spherically-symmetric models, that the collapse 
does not halt once the DM heating rate exceeds the 
H2 cooling rate, and in our present study we confirm 
their results. The key factor that allows collapse to con- 
tinue is the collisional dissociation of H2. This acts as 
a thermostat, preventing the gas temperature from in- 
creasing significantly above 2000 K until all of the H2 is 
destroyed. A simple estimate of the timescale on which 
the H2 is destroyed shows that even in the most extreme 
case, it is comparable to the free-fall time of the gas, 
meaning that the gas can reach much higher densities 
(and hence smaller l ength scales) than anticipated in the 
ISpolvar et al.l (|2008l ) study. As a result, we find no ev- 
idence for the formation of dark stars within our simu- 
lations, although we note that we cannot rule out the 
formation of such objects if they form above densities of 
10 16 cm" 3 , have initial masses of less than 0.01 M Q and 
radii of less than an AU, as this is below our resolution 
limit. 

6.2. Stellar Multiplicity 

Perhaps the most striking finding of our work is the 
large reduction of the level of fragmentation within the 
protostellar accretion disc. Of our four simulations with 
DM, secondary fragmentation occurred in only two cases 
and then only at large radii. Much of the recent focus of 
studies of Population III star formation has been on the 
fragmentation of protos t ellar discs (see e.g . iClark et"aL 
2011b; iGreif et al.l [20TTI iSmith et al.l[201lh . iGreif et al. 
(2012) have shown that fragmentation occurs on scales 
of 10 AU or smaller, and that mergers between proto- 
stars should be common. The possibility of mergers and 
interactions betwee n protostars was also highlighted by 
ISmith et al.l (j2012), who show that the high accretion 
rates experienced by the young protostars lead to large, 
extended, "fluffy" objects that have a higher probability 
of interacting than more conventional protostars. 

Such scenarios are far less likely when the effects of 
DMA are taken into account. In this picture either a 
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Fig. 10. — Column density projections of a slice through the midplane of the disc in the central regions of run H2. The slice is 1000 AU 
thick. In this case a secondary protostar is formed, but at a very large separation from the central protostar. In the second panel, the 
image is centred on the original position of the first protostar. It is currently unclear how the movement of the protostar after its formation 
will affect the DM profile. 



single massive star forms at the centre of the halo, or 
a secondary forms in a wide binary with a separation 
of around 1, 000 AU. Fragmen tation on such scales was 
also seen bv lStacv et al.l (|2010D in lower resolution stud- 
ies which focussed on scales larger than the inner disc 
around the protostar. In our model, when secondary 
fragmentation begins the mass of the pri mary protostar 
is alre ady greater than 10 solar masses. McKee fc Tanl 
(2008) find that for an accretion rate of 10~ 3 M Q yr -1 an 
ionised HII region will form around the protostar after 
it obtains a mass of around 20 M Q , depending upon 
its rotation. Such an ionised region is likely to form 
even earlier in our model due to the additional ionisa- 
tion of hydrogen from DMA. Thus it is possible that at 
the point when secondary stars form the original pro- 
tostar will already be surrounded by an HII region and 
any further growth of the primary will b e determined by 
the ba lance of radiative transfer effects. [Hosokaw a et al.l 
(2012) find that UV radiation from primordial proto- 
stars would photo-evaporate their accretion discs once 
they achieved a mass of around 43 M Q . We note how- 
ever that s imilar claims for pre s ent d ay high-mass star 
formation (jYorke fe Sonnhalterl I2002T) have not borne 
out in detailed 3D calculations (IKrumholz "eFaT][200l 
IPeters et alJl2010HKuiper et al.|[2012l ). 

6.3. Accuracy of DM model 

The biggest assumption in our model is the adoption 
of an analytic DM halo rather than a live halo. No cal- 
culation has yet been able to fully follow the contraction 
and fragmentation of baryonic matter down to AU scales 
with a compara ble resolution in DM. Simulations by 
lAbel et ail (|2002f ) followed the collapse of a minihalo with 
a mass resolution of 1.1 M Q for the DM component and 
confirmed that the DM had a peaked profile to radii as 



small as ~ 1, 500 AU. Previous w ork bv lRipamonti et akl 
(2010) a ndlSpolvar et ail (120 08) adopt the method pro- 
posed bv lBlumenthal et al.l (|1986f ) to model the contrac- 
tion of a DM halo due to baryonic collapse, and upon 
these results we base ou r DM profile. In the Appendix of 
IRipamonti et al.l (|2010D . the applicability of this method 
is discussed and found to be in better agreement with 
the current num eric al find ings than alternative models 
iSteigman et al.l (e.g. I1978D . We are therefore confident 
that our model for increasing the DM density with in- 
creasing baryon density is reasonable. If anything such 
a profile may give a slight overestimate of DM density 
at the centre of the halo. This would further reduce the 
effective of DMA heating, and make it even less likely 
that any "dark stars" are formed. 

A bigger potential uncertainty in our adopted DM pro- 
file, however, is the assumption of spherical symmetry. 
Before the formation of the first sink particle this is a 
reasonable simplification, as the baryon profile itself is 
smooth and centrally concentrated. However, after the 
first sink particle is formed, this rapidly changes within 
the central regions of the halo. A disc quickly devel- 
ops, and strong spiral arm features develop within the 
disc. It is possible that interactions between these high 
asymmetric features and inhomogeneities in the DM den- 
sity distribution could lead to a transfer of energy to the 
DM component and a flattening of the central DM dis- 
tribution, in a similar manner to that thought to occur 
from a rotating bar in disc galaxies (jWeinberg fc Kat"3 
2002). To fully resolve these issues, calculations with a 
live three-dimensional halo component would be needed. 

Another important issue is whether the central proto- 
star remains within the centre of the DM distribution, as 
this affects the amount of DM that can be captured by 
the protostar by particle scattering, which strongly in- 
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fluen c es the lifetime of any dark star phase ()Iocco et al.l 
[20081 iTaoso et al.ll2008at lYoon et alj|200l . In our cal- 
culations, we find that the central protostar moves by 
a significant amount after its formation, as a result of 
the fact that accretion onto the protostar does not occur 
in a perfectly symmetrical fashion, and also because it 
is being acted on by torques from the asymmetrical gas 
distribution in the protostellar accretion disc. We have 
assumed that the DM responds in a similar fashion, and 
that the peak in the DM distribution follows the central 
protostar. However, if an offset were to develop between 
the DM cusp and the baryon peak, dynamical heating of 
the DM would occur and the resulting reduction in the 
DM density would reduce the amount of an nihilation. 

Thi s issue has been addressed recently by iStacv et"aLl 
(|2012j), who performed a calculation in which they re- 
simulated the inner regions of a halo with a live DM halo 
of similar resolution to that of the gas after the formation 
of an initial sink particle. The effects of DMA were not, 
however, included. They formed additional sink particles 
at around 1,000 AU and found that the motion of the 
star-disc system became displaced from the DM density 
peak and that after 5,000 yr there was insufficient DM to 
influence the formation and evolution of the protostars. 
However. lStacv et al.l ((2012) do not resolve protostar for- 
mation in the inner disc around the central object and 
do not include DMA and so cannot comment on frag- 
mentation in this regime. In our runs without DMA, 
secondary fragmentation occurred on timescales of only 
a few hundred years. It is therefore likely that there will 
be sufficient DM during the period in which the inner 
disc is prone to fragmentation to maintain the findings 
of this paper. 

7. CONCLUSIONS 

We have for the first time included the effects of 
dark matter annihilation into 3D simulations of the 
collapse and fragmentation of cosmological minihalos. 
We used the SPH code gadget 2, which has been 
modified to include a time-dependent chemical network 
which includes the effects of dark matter annihilations. 
The dark matter distribution was modelled analyti- 
cally based on the pred ictions of adiabat ic contraction 
( Blumenthal et al.lll986l) and the re sults of lSpolvar et al.l 
( 20081 and lRipamonti et~aTl (|2010t l. We adopted a fidu- 
cial dark matter particle mass of 100 GeV but also ran 
simulations with particles masses an order of magnitude 
larger and smaller to account for uncertainties in our 
knowledge of the true dark matter particle mass and an- 
nihilation cross-section. 



Our main findings are as follows: 

1. Dark matter annihilation does not halt the grav- 
itational collapse of the gas on the scales con- 
sidered here, c ontra ry to the suggestion made by 
ISnolvar et all (l2008h . The reason for this is that 
collisional dissociation of H2 provides an effective 
way to dissipate large amounts of energy even in 
the regime where H2 line cooling is ineffective, and 
hence prevents the gas temperature from increas- 
ing much above 2000 K. The timescale for H2 de- 
struction is typically longer than the free-fall time 
of the gas, even at high gas densities, and so the 
central protostar should be able to reach very high 
densities. This is true even when a dark matter 
particle mass of only 10 GeV is adopted, which 
gives a maximal effect for dark matter annihila- 
tion. It is therefore plausible that a conventional 
star will form, rather than a "dark star" , although 
confirmation of this point awaits a simulation able 
to follow the collapse of the gas all the way up to 
protostellar densities. 

2. Dark matter annihilation prevents subsequent frag- 
mentation in the disc around the central protostar 
within a radius 1,000 AU, regardless of whether 
this central object is a "dark star" or a normal 
protostar. This is primarily due to the dark mat- 
ter annihilation heating raising the temperature of 
the dense gas to a level at which the protostellar ac- 
cretion disc becomes stable. Over time, dark mat- 
ter annihilation will destroy all of the H2 within 
the vicinity of the central protostar, making fur- 
ther fragmentation impossible. At distances larger 
than around 1,000 AU, however, the effects of dark 
matter annihilation are less pronounced and frag- 
mentation can occur in some cases, typically result- 
ing in the formation of a wide binary system. 
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